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In complex systems such as spin systems and protein systems, conventional simulations in 
the canonical ensemble will get trapped in states of energy local minima. We employ the 
generalized-ensemble algorithms in order to overcome this multiple-minima problem. Two well- 
known generalized-ensemble algorithms, namely, muhicanonical algorithm and replica-exchange 
. method, are described. We then present four new generalized-ensemble algorithms as further ex- 

O ' tensions of the two methods. Effectiveness of the new methods are illustrated with a Potts model, 

^ I Lennard- Jones fluid system, and protein system. 

; I. INTRODUCTION 

C/5 . 

In complex systems such as spin systems and protein systems, conventional simulations in the canonical ensemble will 
get trapped in states of energy local minima at low temperatures. We employ the generalized-ensemble algorithms in 
order to overcome this multiple-minima problem (for reviews, see Refs. 3]). In a generalized-ensemble simulation, 
each state is weighted by a non-Boltzmann probability weight factor so that a random walk in potential energy space 
may be realized. The random walk allows the simulation to escape from any energy barrier and to sample much wider 
Q configurational space than by conventional methods. Monitoring the energy in a single simulation run, one can obtain 
^ not only the global-minimum-energy state but also canonical ensemble averages as functions of temperature by the 
single-histogram Q and multiple-histogram |^ reweighting techniques. 
I One of the most well-known generalized-ensemble methods is perhaps multicanonical algoritfim (MUCA) |0] . (The 

^ method is also referred to as entropic sampling 7] , adaptive umbrella sampling of ttie potential energy Q , random 
^\ walk algorith m llOl . and density of states Monte Carlo llj.) MUCA was first introduced to the molecular simulation 
T-H field in Ref. |l2l |. Since then MUCA has been extensively used in many applications in protein and related systems 
(for a review, see, e.g., Ref. ^1,]). 

The replica-exchange method (REM) is another widely used generalized-ensemble algorithm. (Closely related 

. methods were independently developed in Refs. REM is also referred to as multiple Markov chain method 

s ' ^3 and parallel tempering 19|. For recent reviews with detailed references about the method, see, e.g., Refs. UM-) 
REM has also been introduced to protein systems |2]|-[2^. 

Both MUCA and REM are already very powerful, but we have also developed several new generalized-ensemble 
g [ algorithms as further extensions of MUCA and/or REM [23|-|3^. 

I ■ In this article, we first describe the two familiar methods: MUCA and REM. We then present some of our new 
' O ' generalized-ensemble algorithms. The effectiveness of these methods is illustrated with a 2-dimensional Potts model, 
^ Lennard-Jones fluid system, and protein system. 
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II. METHODS 



^ ' In the regular canonical ensemble with a given inverse temperature (3 = l/ksT [ks is the Boltzmann constant), 
5^ , the probability distribution of potential energy E is given by 

Pb{E;T) oi n{E) Wb{E-T) = n{E) e-P^ , (1) 

where n{E) is the density of states. Since the density of states n{E) is a rapidly increasing function of E and the 
Boltzmann factor Wb{,E]T) decreases exponentially with E, the probability distribution Pb{E;T) has a bell- like 
shape in general. A Monte Carlo (MC) simulation based on the Metropolis algorithm generates states in the 
canonical ensemble with the following transition probability from a state x with energy iJ to a state x' with energy 
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E': 

wi. ^ x') = min (l, ^^^) - min (l, e-^^^'-^^) . (2) 

However, it is very difficult to obtain canonical distributions at low temperatures with this conventional Metropolis 
algorithm. This is because the thermal fluctuations at low temperatures are small and the simulation will certainly 
get trapped in states of energy local minima. 

In the "multicanonical ensemble" on the other hand, the probability distribution of potential energy is defined 
as follows so that a uniform flat distribution of E may be obtained: 

Pmu{E) cx n{E) W,nu{E) = constant . (3) 

Hence, the multicanonical weight factor Wmu{E) is inversely proportional to the density of states, and the Metropolis 
criterion for the multicanonical MC simulations is based on the following transition probability: 

Because the MUCA weight factor Wmu{E) is not a priori known, however, one has to determine it for each system 
by iterations of trial simulations. 

After the optimal MUCA weight factor is obtained, one performs a long MUCA simulation once. By monitoring the 
potential energy throughout the simulation, one can find the global-minimum-energy state. Moreover, by using the 
obtained histogram N^a{E) of the potential energy distribution Pmu(-E'), the expectation value of a physical quantity 
A at any temperature T = \/k^[3 can be calculated from 

A{E) n{E) e-f^^ 

<A>T = Af^ , (5) 



e -'^^ 



where the best estimate of the density of states is given by the single-histogram reweighting techniques (see Eq. Q) 
i: 

The system for replica- exchange method (REM) 0, consists of M non-interacting copies, or replicas, of the 
original system in canonical ensemble at M different temperatures Tm (m = 1, • • • , M). We arrange the replicas so 
that there is always one replica at each temperature. Then there is a one-to-one correspondence between replicas 
and temperatures. Let X — | - • • ,Xm, ' ' ' } ^^^'^^ ^ state in this generalized ensemble. Here, the superscript 

i and the subscript m in Xm label the replica and the temperature, respectively. A simulation of REM is then 
realized by alternately performing the following two steps. Step 1: Each replica in the canonical ensemble at a 
fixed temperature is simulated simultaneously and independently for a certain number of MC steps. Step 2: A 
pair of replicas, say i and j, which are at neighboring temperatures, say Tm and Tm+i, respectively, are exchanged: 

X ~ I • • • ,Xm, - ■ • , x^+i, • • • I ^ X' = I • • • , Xm , • • • , Xm+1 ''''}■ ^hc transition probability of this replica exchange 
is given by the following Metropolis criterion: 

w{X ^ X') =iam{l,e-^) , (7) 

where 

A^(/3™+i-A„)(i?(gW) . (8) 

From the results of a long REM production run, one can obtain the canonical ensemble average of a physical 
quantity A as a function of temperature from Eq. ||SJ|, where the density of states is given by the multiple- histogram 
reweighting techniques 5] as follows. Let Nm{E) and Um be respectively the potential-energy histogram and the 
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total number of samples obtained at temperature T™ = l/fce/Jm (™ = 1, • ■ • , M). The best estimate of the density 
of states is then given by |^ 



M 



J2 97n Nra{E) 

^(E) = -tP^ , (9) 



M 

97n nm e- 

m—1 



where 



^ n{E) e-P-^ . (10) 



Here, ffm = 1 + 2T„i, and Tm is the integrated autocorrelation time at temperature T^. Note that Eqs. Q and l|10() 
are solved self-consistently by iteration Q to obtain the dimensionless Helmholtz free energy and the density of 
states n{E). 

We now introduce new generalized-ensemble algorithms that combine the merits of MUCA and REM. In the replica- 
exchange multicanonical algorithm (REMUCA) [2^f3lj | we first perform a short REM simulation (with M replicas) to 
determine the MUCA weight factor and then perform with this weight factor a regular MUCA simulation with high 
statistics. The first step is accomplished by the multiple- histogram reweighting techniques j^]. Let Nm{E) and n„i be 
respectively the potential-energy histogram and the total number of samples obtained at temperature Tm — l/fce/Sm 
of the REM run. The density of states n{E), or the inverse of the MUCA weight factor, is then given by solving 
Eqs. lini and H1U|I self-consistently by iteration |^. The formulation of REMUCA is simple and straightforward, but 
the numerical improvement is great, because the weight factor determination for MUCA becomes very difficult by the 
usual iterative processes for complex systems. 

While multicanonical simulations are usually based on local updates, a replica-exchange process can be considered 
to be a global update, and global updates enhance the sampling further. Here, we present a further modification 
of REMUCA and refer to the new method as multicanonical replica- exchange method (MUCAREM) (29.. .31.j . In 
MUCAREM the final production run is not a regular multicanonical simulation but a replica-exchange simulation 
with a few replicas in the multicanonical ensemble. Because multicanonical simulations cover much wider energy 
ranges than regular canonical simulations, the number of required replicas for the production run of MUCAREM 
is much less than that for the regular REM, and we can keep the merits of REMUCA (and improve the sampling 
further). The details of REMUCA and MUCAREM can be found in Ref. 0] 

Besides canonical ensemble, MC simulations in isobaric-isothcrmal ensemble are also extensively used. This 
is because most experiments are carried out under the constant pressure and constant temperature conditions. The 
distribution Pnpt{E, V) for E and V is given by 

PnpAE, V) = n{E, V)e-'^''" . (11) 

Here, the density of states n{E, V) is given as a function of both E and V, and H is the "enthalpy": 

H = E + PoV, (12) 

where Pq is the pressure at which simulations are performed. This ensemble has bell-shaped distributions in both E 
and V. 

We now introduce the idea of the multicanonical technique into the isobaric-isothermal ensemble MC method 
and refer to this generalized-ensemble algorithm as the multibaric-multithermal algorithm |32j| . This MC simulation 
performs random walks in volume space as well as in potential energy space. 

In the multibaric-multithermal ensemble, each state is sampled by a weight factor Wi-aht{E,V) = 
exp{— /3o-ffmbt(^', V)} {H-caht is referred to as the multibaric-multithermal enthalpy) so that a uniform distribution in 
both potential energy and volume is obtained: 



mbt 



{E, V) = n{E, U)MVbt {E, V) = constant . (13) 



We call Wn±,t{E,V) the multibaric-multithermal weight factor. 

In order to perform the multibaric-multithermal MC simulation, we follow the conventional isobaric-isothermal MC 
techniques [33 ■ In this method, we perform Metropolis sampling on the scaled coordinates Si = L^^Vi (r^ are the 
real coordinates) and the volume V (here, the particles are placed in a cubic box of a side of size L = -v^). The 
trial moves of the scaled coordinates from Si to s'i and of the volume from V to V are generated by uniform random 
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FIG. 1: (a) Probability distributions of energy of 2-dimensional 10-state Potts model at three temperatures; T = 0.6000, 
0.7026, and 0.8000, and (b) average energy and (c) specific heat as functions of inverse temperature 13 — 1/T. The results 
were obtained from a multicanonical MC simulation. For (b) and (c), the results from three methods of multicanonical weight 
factor determination are superimposed. Berg stands for Berg's method [39|. W-L for Wang-Landau's method ,10,| and MR for 
MUCAREM H 



numbers. The enthalpy is accordingly changed from H{E{s^^\ V), V) to H'{E{s''^^\ V'), V) by these trial moves. 
The trial moves will be accepted with the probability 

=min(l,exp[-/3o{-ff'-i/-iVfcsToln(l/7F)}]) , (14) 

where N is the total number of particles in the system. 

Replacing H by iJmbtj we can perform the multibaric-multithermal MC simulation. The trial moves of Si and V 
are generated in the same way as in the isobaric-isothermal MC simulation. The multibaric-multithermal enthalpy 
is changed from H,^ht{E{s'-^\V),V) to H^-^^^{E{s"^'^\V'),V') by these trial moves. The trial moves will now be 
accepted with the probability 

wix ^ x') = min (1, exp[-/3o{i/Ut - ^mbt - NkeTo HV /V)}]) . (15) 

While MUCA yields a flat distribution in potential energy and performs a random walk in potential energy space, 
we can, in principle, choose any other variable and induce a random walk in that variable. One such example is the 
multi-overlap algorithm |33j| . Here, we choose a protein system and define the overlap in the space of dihedral angles 
by, as it was already used in [36l| . 

q = (n — d)/n , (16) 
where n is the number of dihedral angles and d is the distance between configurations defined by 

n 

d = \\v-v'\\ = -y^da{v^,vl) . (17) 

1=1 

Here, Vi is our generic notation for the dihedral angle i, — tt < < tt, and is the vector of dihedral angles of the 
reference configuration. The distance da{vi,v[) between two angles is defined by 

da{vi,vl) = min{\vi - Vi\,2Tr - \vi - v'i\) . (18) 

We want to simulate the system with weight factors that lead to a flat distribution in the dihedral distance d, and 
hence to a random walk process in d: 

d < dmin d > dmax and back . (19) 

Here, dmin is chosen sufficiently small so that one can claim that the reference configuration has been reached. The 
value of dmax has to be sufficiently large to introduce a considerable amount of disorder. 

Moreover, we can define a weight factor that leads to a random-walk process between two configurations. This 
multi-overlap simulation allows a detailed study of the transition states between the two configurations, whereas a 
random walk in energy space of a regular MUCA simulation may miss the transition state (see Ref. for details). 
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FIG. 2; Probability distributions of energy for the 2-dimensional 10-state Potts model; (a) the results of REM simulation with 
32 replicas and (b) and (c) iterations of MUCAREM simulations with 8 replicas. 
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FIG. 3: Probability distribution of energy during the iterative process of multicanonical weight factor determination for the 
2-dimensional 10-state Potts model. The results after 300^00 MC sweeps, 500,000 MC sweeps, and 1,000,000 MC sweeps are 
superimposed, (a) MUCAREM |2|], (b) Berg's method |3|], and (c) Wang-Landau's method [lo(|. 



III. RESULTS 



We now present the results of our simulations based on the algorithms described in the previous section. 

The first example is a spin system. We studied the 2-dimensional 10-state Potts model [s^- The lattice size was 
34 X 34. This system exhibits a first-order phase transition _38j. In Fig. ^ we show the probability distributions 
of energy at three tempeartures (above the critical temperature Tc, at Tq, and below Tc) and average energy and 
specific heat as functions of inverse temperature. All these results imply that the system indeed undergoes a first-order 
phase transition. 

Iterations of MUCAREM (and REMUCA) can be used to obtain an optimal MUCA weight factor In Fig.|3 
we show the results of our MUCA weight factor determination by MUCAREM. We first made a REM simulation 
of 10,000 MC sweeps (for each rephca) with 32 replicas (Fig. I2a)). Using the obtained energy distributions, we 
determined the (preliminary) MUCA weight factor, or the density of states, by the multiple-histogram reweighting 
techniques of Eqs. 10 and (|10|l . Because the trials of replica exchange are not accepted near the critical temperature 
for first-order phase transitions, the probability distributions in Fig.[2Ia) for the energy range from ^ —1.5 to ^ —1.0 
fails to have sufficient overlap, which is required for successful application of REM. This means that the MUCA weight 
factor in this energy range thus determined is of "poor quality." With this MUCA weight factor, however, we made 
iterations of three MUCAREM simulations of 10,000 MC sweeps (for each replica) with 8 replicas (Fig. |2Ib) for the 
first iteration and Fig.[2Ic) for the third iteration). In Fig. |2Ib) we see that the distributions are not completely flat, 
reflecting the poor quality in the phase-transition region. This problem is rapidly rectifled as iterations continue, and 
the distributions are completely flat in Fig. |2Ic) , which gives an optimal MUCA weight factor in the entire energy 
range by the multiple-histogram reweighting techniques. 

Besides MUCAREM the methods of Berg and Wang-Landau |0 are also effective for the determination of the 
MUCA weight factor (or the density of states). In Fig. Owe compare how fast these three methods converge to yield 
an optimal MUCA weight factor, or a flat distribution. Each flgure shows three curves superimposed that correspond 
to the (hot-start) MUCA simulation with the weight factor that was "frozen" after 300,000 MC seeps, 500,000 MC 
sweeps, and 1,000,000 MC sweeps of iterations of the weight factor determination. While the results of MUCAREM 
and Berg's method are similar, those of Wang-Landau method have quite different behavior. After 300,000 MC 
sweeps, MUCAREM and Berg's method gives reasonably flat distribution for E > —1.0 and E > —1.3, respectively, 
whereas Wang-Landau method gives a distribution that is quite spiky (but already covers low-energy regions) . After 
500,000 MC sweeps, MUCAREM essentially gives a flat distribution in the entire energy range and Berg's method 
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FIG. 4: (a) The probability distribution Pnpt{E* /N,V* /N) in the isobaric-isothermal simulation at (T*,P*) = {To,Po) 
(2.0,3.0) and (b) the probability distribution P^nhtiE* /N,V* /N) in the multibaric-multithermal simulation. 
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FIG. 5: The time series of V* /N from (a) the conventional isobaric-isothermal MC simulations at (T*, P*) — (2.0, 2.2) and at 
(T*,P*) = (2.0,3.8) and (b) the multibaric-multithermal MC simulation. 



for E > —1.7, while Wang-Landau method also covers the entire range (though still very rugged). After 1,000,000 
MC sweeps, the three methods all give reasonably flat distributions. Details will be published elsewhere 

We now present the results of our multibaric-multithermal simulation. We considered a Lennard- Jones 12-6 potential 
system. We used 500 particles {N = 500) in a cubic unit cell with periodic boundary conditions. The length and the 
energy are scaled in units of the Lennard-Jones diameter a and the minimum value of the potential e, respectively. 
We use an asterisk (*) for the reduced quantities such as the reduced length r* — r/a, the reduced temperature 
T* — k-QT/e, the reduced pressure P* = Pcr'^/e, and the reduced number density p* ~ pa^ [p = N/V). 

We started the iterations of the multibaric-multithermal weight factor determination from a regular isobaric- 
isothermal simulation at Tq = 2.0 and Pq = 3.0. In one MC sweep we made the trial moves of all particle coordinates 
and the volume {N + 1 trial moves altogether). For each trial move the Metropolis evaluation of Eq. (|15|l was made. 
Each iteration of the weight factor determination consisted of 100,000 MC sweeps. In the present case, it was required 
to make 12 iterations to get an optimal weight factor Wnibt(P, V). We then performed a long multibaric-multithermal 
MC simulaton of 400,000 MC sweeps with this Wmbt(£^, V). 

Figure 01 shows the probability distributions of E* /N and V* /N. Figure Efa) is the probability distribution 
P npt{E* /N,V* /N) from the isobaric-isothermal simulation first carried out in the process (i.e., Tg* = 2.0 and 
Pq — 3.0). It is a bell-shaped distribution. On the other hand, Fig. E^b) is the probability distribution 
Pniht{E* /N,V* /N) from the multibaric-multithermal simulation finally performed. It shows a flat distribution, and 
the multibaric-multithermal MC simulation indeed sampled the configurational space in wider ranges of energy and 
volume than the conventional isobaric-isothermal MC simulation. 

Figure |S1 shows the time series of V* /N. In Fig. EJa) we show the results of the conventional isobaric-isothermal 
simulations at {T*,P*) — (2.0,2.2) and (2.0,3.8), while in Figure|SIb) we give those of the multibaric-multithermal 
simulation. The volume fluctuations in the conventional isobaric-isothermal MC simulations are only in the range 
of V*/N = 1.3 1.4 and V*/N = 1.5 ~ 1.6 at P* = 3.8 and at P* = 2.2, respectively. On the other hand, the 
multibaric-multithermal MC simulation performs a random walk that covers even a wider volume range. 
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FIG. 6: (a) Average potential energy per particle < E* /N >npt and (b) average density < p* >npt at various temperature 
and pressure values. Filled circles: Multibaric-multithermal MC simulations. Open squares: Conventional isobaric-isothermal 
MC simulations. Solid line: Equation of states calculated by Johnson et al. [4(1 |. Broken line: Equation of states calculated by 
Sun and Teja 




FIG. 7: (a) Reference configuration 1 and (b) reference configuration 2. Only backbone structures are shown. The N-terminus 
is on the left-hand side and the C-terminus on the right-hand side. The dotted lines stand for hydrogen bonds. The figures 
were created with RasMol |4^. 

We calculated the ensemble averages of potential energy per particle, < E*/N >npt, and density, < p* >npt, at 
various temperature and pressure values by the reweighting techniques. They are shown in Fig. |S| The agreement 
between the multibaric-multithermal data and isobaric-isothermal data are excellent in both < E*/N >npt and 

< P* >NPT- 

The important point is that we can obtain any desired isobaric-isothermal distribution in wide temperature and 
pressure ranges (T* — 1.6 ~ 2.6, P* — 2.2 ~ 3.8) from a single simulation run by the multibaric-multithermal 
MC algorithm. This is an outstanding advantage over the conventional isobaric-isothermal MC algorithm, in which 
simulations have to be carried out separately at each temperature and pressure, because the reweighting techniques 
based on the isobaric-isothermal simulations can give correct results only for narrow ranges of temperature and 
pressure values. 

The third example is a system of a biopolymer. A brain peptide Met-enkephalin has the amino-acid sequence 
Tyr-Gly-Gly-Phe-Met. We fix the peptide-bond dihedral angles oj to 180°, which implies that the total number of 
variable dihedral angles is n = 19. We neglect the solvent effects as in previous works. The low-energy configurations 
of Met-enkephalin in the gas phase have been classified into several groups of similar structures [i^. l4a| . Two 
reference configurations, called configuration 1 and configuration 2, are used in the following and depicted in Fig.EI 
Configuration 1 has a /3-turn structure with hydrogen bonds between Gly-2 and Met-5, and configuration 2 a /3-turn 
with a hydrogen bond between Tyr-1 and Phe-4 :43j. Configuration 1 corresponds to the global- minimum-energy 
state and configuration 2 to the second lowest-energy state. The distance between the two configurations is o? = 6.62 
and the values of the potential energy (ECEPP/2 44]) for configuration 1 and configuration 2 are —10.72 kcal/mol 
and —8.42 kcal/mol, respectively. 

We analyze the free-energy landscape |43i from the results of our multi-overlap simulation at 300 K that performs 
a random walk between configurations 1 and 2. We study the landscape with respect to some reaction coordinates 
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FIG. 8: (a) Free-energy landscape of Met-enkephalin a,t T — 250 K with respect to rms distances (A) from the two reference 
configurations, F{ri, r2). The labels Ai and Bi indicate the positions for the local-minimum states at T = 250 K that originate 
from the reference configuration 1 and the reference configuration 2, respectively. The label C stands for the saddle point that 
corresponds to the transition state, (b) The transition state, C, between reference configurations 1 and 2. See the caption of 
Fig.[7|for details. 



TABLE I: Free energy, internal energy, entropy multiplied by temperature at T = 250 K (all in kcal/mol) at the two local- 
minimum states (Ai and Bi) and the transition state (C) in Fig. |HJa). The rms distances, ri and r2, are in A. 



Coordinate (ri,r2) 


F 


U 


~TS 


Ai (1.23, 4.83) 





-5.4 


5.4 


Bi (4.17, 2.43) 


1.0 


-3.5 


4.5 


C (3.09, 4.05) 


2.2 


-0.8 


3.0 



(and hence it should be called the potential of mean force). In order to study the transition states between reference 
configurations 1 and 2, we first plotted the free-energy landscape with respect to the dihedral distances di and d2 of 
Eq. (|17|l . However, we did not observe any transition saddle point. A satisfactory analysis of the saddle point becomes 
possible when the root-mean-square (rms) distance (instead of the dihedral distance) is used. Fig ure |H1 shows contour 
lines of the free energy reweighted to T = 250 K, which is close to the folding temperature ([IJSI)- Here, the free 
energy F{ri,r2) is defined by 

F(ri,r2) = ~kBT\nP{ri,r2) , (20) 

where ri and r2 are the rms distances from the reference configuration 1 and the reference configuration 2, respectively, 
and P{ri,r2) is the (reweighted) probability at T = 250 K to find the peptide with values ri,r2. The probability 
was calculated from the two-dimensional histogram of bin size 0.06 Ax 0.06 A. The contour lines were plotted every 
2kBT (= 0.99 kcal/mol for T = 250 K). 

Note that the reference configurations 1 and 2, which are respectively located at (ri,r2) = (0,4.95) and (4.95,0), 
are not local minima in free energy at the finite temperature (T = 250 K) because of the entropy contributions. The 
corresponding local-minimum states at Ai and Bi still have the characteristics of the reference configurations in that 
they have backbone hydrogen bonds between Gly-2 and Met-5 and between Tyr-1 and Phe-4, respectively. 

The transition state C in Fig.|S{a) should have intermediate structure between configurations 1 and 2. In Fig.lS^b) 
we show a typical backbone structure of this transition state. We see the backbone hydrogen bond between Gly-2 
and Phe-4. This is precisely the expected intermediate structure between configurations 1 and 2, because going from 
configuration 1 to configuration 2 we can follow the backbone hydrogen-bond rearrangements: The hydrogen bond 
between Gly-2 and Met-5 of configuration 1 is broken, Gly-2 forms a hydrogen bond with Phe-4 (the transition state), 
this new hydrogen bond is broken, and finally Phe-4 forms a hydrogen bond with Tyr-1 (configuration 2). 

In Ref. m the low -energy conformations of Met-enkephalin were studied in detail and they were classified into 
several groups of similar structures based on the pattern of backbone hydorgen bonds. It was found there that below 
T = 300 K there are two dominant groups, which correspond to configurations 1 and 2 in the present article. Although 
much less conspicuous, the third most populated structure is indeed the group that is identified to be the transition 
state in the present work. 

In tabled we list the numerical values of the free energy, internal energy, and entropy multiplied by temperature at 
the two local-minimum states (Ai and Bi in Fig. (SJa)) and the transition state (C in Fig.|S{a)). Here, the internal 
energy U is defined by the (reweighted) average ECEPP/2 44] potential energy: 



f/(ri,r2) =< E{ri,r2) >t 



(21) 
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The entropy S was then calculated by 

S{n,r2) ^ ^[U(n,r2) - F{rur2)] • (22) 

The free energy was normalized so that the value at Ai is zero. 

The state Ai can be considered to be "deformed" configuration f, and Bi deformed configuration 2 due to the 
entropy effects, whereas C is the transition state between Ai and Bi. Among these three points, the free energy F 
and the internal energy U are the lowest at Ai, while the entropy contribution —TS is the lowest at C. The free 
energy difference AF, internal energy difference AU, and entropy contribution difference —TAS are f .0 kcal/mol, f .9 
kcal/mol, and —0.9 kcal/mol between Bi and Ai, 2.2 kcal/mol, 4.6 kcal/mol, and —2.4 kcal/mol between C and Ai, 
and f .2 kcal/mol, 2.7 kcal/mol, and — f.5 kcal/mol between C and Bi. Hence, the internal energy contribution and 
the entropy contribution to free energy are opposite in sign and the magnitude of the former is roughly twice as that 
of the latter at this temperature. 



IV. CONCLUSIONS 



In this article we have described the formulations of the two well-known generalized-ensemble algorithms, namely, 
multicanonical algorithm (MUCA) and replica-exchange method (REM). We then introduced four new generalized- 
ensemble algorithms as further extensions of the above two methods, which we refer to as replica-exchange mul- 
ticanonical algorithm (REMUCA), multicanonical replica-exchange method (MUCAREM), multibaric-multithcrmal 
algorithm, and multi-overlap algorithm. 

With these new methods available, we believe that we now have working simulation algorithms for spin systems 
and biomolecular systems. 
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